/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		segment.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date:  2020-01-17 16:54
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:23
*
*   Description:
*
=============================================================================*/
#ifndef I_SEGMENT_H
#define I_SEGMENT_H

#include "segment_shape_function.hpp"
#include "shape_interface.hpp"

namespace geoxfem {
namespace geometry {

template <typename DimensionType, typename degree>
struct Segment : ShapeInterface<DimensionType, degree> {
  DEFINE_KEY(Segment)
  using type = boost::geometry::linestring_tag;

  using Shape      = ShapeInterface<DimensionType, degree>;
  using point_type = typename traits::point_type<DimensionType>::type;
  using Gauss      = model::GaussPoint<DimensionType>;

  static constexpr size_t dimension = traits::dimension<DimensionType>::value;
  /*constructors and destructor{{{*/

  /**
   * @brief constructor
   *
   * @param _listPoint
   */
  Segment(const std::vector<point_type>& listPoint) : Shape{listPoint} {}
  /* Segment(const std::vector<point_type>&& listPoint) : Shape{listPoint} {} */
  Segment()          = default;
  virtual ~Segment() = default;
  /*}}}*/
  /*misc{{{*/

  Segment<DimensionType, degree>& toSegment() final override {
    return *this;
  }
  const Segment<DimensionType, degree>& toSegment() const final override {
    return *this;
  }
  /**
   * @brief getNumberofPoint
   *
   * @return
   */
  size_t getNumberOfPoint() const final override {
    return degree::value + 1;
  }
  /**
   * @brief get Class Name
   *
   * @return
   */
  std::string getClassName() const final override {
    return std::string("Segment") + Shape::poly_degree[degree::value]
           + std::to_string(dimension) + "D";
  }
  /**
   * @brief get Topology Dimension of current shape
   *
   * @return
   */
  size_t getTopologyDimension() const final override {
    return 1;
  }
  /**
   * @brief getVTKTypeID
   *
   * @return
   */
  size_t getVTKTypeID() const final override {
    return getVTKTypeIDImpl<degree>();
  }
  template <typename d>
  size_t getVTKTypeIDImpl() const {
    return 51;
  }
  template <>
  size_t getVTKTypeIDImpl<mpl1>() const {
    return 3;
  }
  template <>
  size_t getVTKTypeIDImpl<mpl2>() const {
    return 21;
  }
  template <>
  size_t getVTKTypeIDImpl<mpl3>() const {
    return 35;
  }

  void plot(io::DataPloter&    plt,
            const std::string& seg_cfg,
            const std::string& gp_cfg) const final override {
    const _real x0 = this->getHeadPoint().x();
    const _real y0 = this->getHeadPoint().y();
    const _real z0 = this->getHeadPoint().z();
    const _real x1 = this->getTailPoint().x();
    const _real y1 = this->getTailPoint().y();
    const _real z1 = this->getTailPoint().z();
    plt.plotLine(x0, y0, z0, x1, y1, z1, seg_cfg);
    if(!gp_cfg.empty())
      for(auto& gp : gpsIn(*this))
        gp.plot(plt, gp_cfg);
    /* for(auto vertex : neighborsIn(*this, 0)) */
    /*   vertex->plot(plt); */
    return;
  }
  /*}}}*/
  /*shape function using reference coordinates{{{*/

  /** calculate standard shape function using local coordinates
   *
   * @param lcoord
   * @return
   */
  VectorXR getStdNByLocal(const VectorXR& lcoord) const final override {
    return segment_shape_function<dimension, degree::value>::getStdNByLocal(
      lcoord, *this);
  }
  /** calculate standard shape function derivation dNdlx using local
   * coordinates
   *
   * @param lcoord
   * @return
   */
  MatrixXR getStddNdlxByLocal(const VectorXR& lcoord) const final override {
    return segment_shape_function<dimension, degree::value>::getStddNdlxByLocal(
      lcoord, *this);
  }
  /*}}}*/
  /*shape function using global coordinates{{{*/

  /** also you should get shape function N with global coordinate gcoord
   *
   * @param gcoord
   * @return
   */
  VectorXR getStdNByGlobal(const VectorXR& gcoord) const final override {
    return segment_shape_function<dimension, degree::value>::getStdNByGlobal(
      gcoord, *this);
  }
  /** then global coordinate can be yield by function above
   *
   * @param gcoord
   * @return
   */
  VectorXR global2local(const VectorXR& gcoord) const final override {
    return segment_shape_function<dimension, degree::value>::global2local(
      gcoord, *this);
  }
  /*}}}*/
  /*gauss point related{{{*/

  /**
   * @brief build standard GaussPoints
   */
  void getStdGaussPoints(std::vector<Gauss>& gp_list) final override {
    size_t   num_gp = model::i_default_gauss_num::set_get().at("segment");
    MatrixXR lcoord(dimension, num_gp);
    VectorXR lweight(num_gp);

    switch(num_gp) {
    case 1:
      lcoord << r0,  //
        r0;
      lweight << r2;
      break;
    case 2:
      lcoord << -r1 / root3, r1 / root3,  //
        r0, r0;
      lweight << r1, r1;
      break;
    case 3:
      lcoord << -root3 / root5, r0, root3 / root5,  //
        r0, r0, r0;
      lweight << static_cast<_real>(5. / 9.),  //
        static_cast<_real>(8. / 9.),           //
        static_cast<_real>(5. / 9.);
      break;
    case 4:
      lcoord << static_cast<_real>(-.861136311594053),
        static_cast<_real>(-.339981043584856),
        static_cast<_real>(.339981043584856),
        static_cast<_real>(.861136311594053),  //
        r0, r0, r0, r0;
      lweight << static_cast<_real>(.347854845137454),
        static_cast<_real>(.652145154862546),
        static_cast<_real>(.652145154862546),
        static_cast<_real>(.347854845137454);
      break;
    case 5:
      lcoord << static_cast<_real>(.906179845938664),
        static_cast<_real>(.538469310105683), r0,
        static_cast<_real>(-.538469310105683),
        static_cast<_real>(-.906179845938664),  //
        r0, r0, r0, r0, r0;
      lweight << static_cast<_real>(.236926885056189),
        static_cast<_real>(.478628670499366),
        static_cast<_real>(.568888888888889),
        static_cast<_real>(.478628670499366),
        static_cast<_real>(.236926885056189);
      break;
    case 6:
      lcoord << static_cast<_real>(-.932469514203152),
        static_cast<_real>(-.661209386466265),
        static_cast<_real>(-.238619186083197),
        static_cast<_real>(.238619186083197),
        static_cast<_real>(.661209386466265),
        static_cast<_real>(.932469514203152),  //
        r0, r0, r0, r0, r0, r0;
      lweight << static_cast<_real>(.171324492379170),
        static_cast<_real>(.360761573048139),
        static_cast<_real>(.467913934572691),
        static_cast<_real>(.467913934572691),
        static_cast<_real>(.360761573048139),
        static_cast<_real>(.171324492379170);
      break;
    default:
      throw "illegal gauss point number!";
      break;
    }
    const _real length = boost::geometry::length(*this);
    for(size_t i = 0; i != num_gp; ++i) {
      VectorXR gcoord  = this->local2global(lcoord.col(i));
      _real    gweight = lweight[i] / r2 * length;
      gp_list.emplace_back(Gauss{gcoord, gweight});
    }
    return;
  }
  /*}}}*/
  /*topology connectivity{{{*/

  std::vector<std::unique_ptr<Shape>> createFacet() const final override {
    std::vector<std::unique_ptr<Shape>> v;
    for(auto& point : pointsIn(*this))
      v.emplace_back(
        Shape::getFactory().create("Vertex", {point.shallowClone()}));

    return v;
  }
  /*}}}*/
};

}  // namespace geometry
}  // namespace geoxfem

GEOXFEM_REGISTER_TEMPLATIZED(geoxfem::geometry::Segment);

#endif /* I_SEGMENT_H */
